rm(list=ls()) 
library(data.table)
library(stargazer)
library(lfe)
library(dplyr)
library(tidyr)
library(haven)

# import temperature: 
tempdf=read.csv("Derived Data/weatherdf.csv",header=TRUE,stringsAsFactors = FALSE)

tempvarsavg="ay_mean_snow  +aymean_high_num_snow+aymean_severe_num_snow+ay_avg_mean_max + ay_avg_mean_min  + ay_avg_bint80  + ay_avg_bint90  + ay_avg_bint100  + ay_avg_bint20 + ay_avg_bint10+ ay_avg_bint0+ ay_mean_prcp +aymean_high_num_prcp+aymean_severe_num_prcp"

districtcovariates="pct_hisp_2005:factor(year) + pct_black_2005:factor(year) + pct_male_2005:factor(year) + pct_white_2005:factor(year) + enroll_2005:factor(year) + pct_ell_2005:factor(year) + pct_speced_2005:factor(year) + pct_FRPM_2005:factor(year) + missing_pct_ell_2005 + missing_pct_speced_2005 + missing_pct_FRPM_2005 "



acs="pctemploy2005:factor(year)  + pctsinglemother2005:factor(year) + pctmanufac2005:factor(year) + pctutility2005:factor(year) + pctbachelorhigher2005:factor(year)"
# import acs

acsdf=read.csv("Derived Data/acs2009_2018.csv",stringsAsFactors = FALSE, header=TRUE)
acsdf$pctemploy=as.numeric(acsdf$pctemploy); acsdf$pctlaborforce=as.numeric(acsdf$pctlaborforce); acsdf$pctutility=as.numeric(acsdf$pctutility); acsdf$pctmanufac=as.numeric(acsdf$pctmanufac)
acsdf$leaid=acsdf$LEAID; acsdf$LEAID=NULL

############### Table A.6 
distance=40


# instrument: 
instname=paste("bartikleave_",distance,sep="")
# weights:
weightval=paste("totalprod_fuel",distance,sep="")




dfwbgraw=fread("Derived Data/df_bwg.csv")

#40km
distance=40
filetoimport=paste("Derived Data/dfbartik2005_",distance,"kmanalysis_annual.csv",sep="")

dfbartik=fread(filetoimport)

dfbartik=merge(dfbartik,tempdf,by=c("leaid","year"),all.x=TRUE)


dfbartik$cs_mn_wbg=NULL



dfwbg=merge(dfbartik,dfwbgraw,by=c("year","grade","subject","leaid"))

dfwbg=merge(dfwbg,acsdf,by=c("leaid","year"),all.x=TRUE)


dfwbgfinal=dfwbg[(dfwbg$totgyb_blk+dfwbg$totgyb_wht)>=50 & dfwbg$totgyb_wht>=25 &dfwbg$totgyb_blk>=25,]


bartik_dist=as.formula(paste("cs_mn_wbg~ lagcs_mn_wht + lagcs_mn_blk  + cs_mn_all_lag_cohort + ", districtcovariates,  " |factor(year) + factor(subject)  + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
resultbartik_dist=felm(bartik_dist,dfwbgfinal,weights=dfwbgfinal$totalprod_fuel40,na.action=na.omit)

bartik_dist_time=as.formula(paste("cs_mn_wbg~ lagcs_mn_wht + lagcs_mn_blk + ", districtcovariates, "+", acs ,"|factor(year)  + factor(grade)+factor(leaid) + factor(subject) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
resultbartik_dist_time=felm(bartik_dist_time,dfwbgfinal,weights=dfwbgfinal$totalprod_fuel40,na.action=na.omit)


bartik_distweather=as.formula(paste("cs_mn_wbg~ lagcs_mn_wht + lagcs_mn_blk+", districtcovariates, "+", acs ,"+", tempvarsavg,"|factor(year)  + factor(grade)+factor(leaid) + factor(subject) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
resultbartik_dist_weather=felm(bartik_distweather,dfwbgfinal,weights=dfwbgfinal$totalprod_fuel40,na.action=na.omit)

stargazer(  resultbartik_dist, resultbartik_dist_time, resultbartik_dist_weather,keep="avgpm25_9",out="output/TableA6.tex",digits=4 )


